################################################################################# 
#                     Load Packages                                            #
################################################################################# 

packages = c("brglm","ggcorrplot","interplot", "stargazer", "tidyverse",
             "here","haven","mice","modelsummary", "sandwich", "lmtest",
             "MASS", "ROCR", "generalhoslem", "brant", "gridExtra",
             "erer", "margins")
#Warning! This will install R packages in your system if not yet installed. 
#load or install & load all packages 
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

################################################################################# 
#                     Prepare Data                                             #
################################################################################# 

#####Inport Data ######
LEDs_data_f <- read.csv(here::here("replication_data_afs.csv")) 

################################################################################# 
#                     Analysis                                             #
################################################################################# 


#Table 1
no_missing <- LEDs_data_f %>% 
  ungroup() %>% dplyr::select(wdldownes, wl,logLER,wldownes) %>%
  drop_na() %>%
  mutate(wdldownes = as.numeric(wdldownes),
         wl = as.numeric(wl),
         logLER = as.numeric(logLER),
         wldownes = as.numeric(wldownes),)

stargazer(cor(no_missing))

#Table 2: Models

m1<- polr(as.factor(wdldownes) ~ relative_corr +
            initally + targally + v2x_libdem+ v2x_ex_military + concap + capasst +
            qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
          data = LEDs_data_f, method = c("probit"))
summary(m1)
ocME(w = m1)
m1sr<- coeftest(m1, 
                vcov = sandwich(x =m1, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m2<- glm(wl ~ relative_corr +
           initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m2)
m2sr<- coeftest(m2, 
                vcov = sandwich(x =m2, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

m3<- glm(wldownes ~ relative_corr +
           initally + targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f,  family = binomial(link = "probit"))
summary(m3)
m3sr<- coeftest(m3, 
                vcov = sandwich(x = m3, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

#with LEDs

m4 <- lm(logLER ~ relative_corr +
           + initally +targally + v2x_libdem+ v2x_ex_military + concap+ capasst+
           qualrat+ terrain+ straterr+ strat1+ strat2+ strat3+ strat4, 
         data = LEDs_data_f) 
summary(m4)
m4sr<- coeftest(m4, 
                vcov = vcovHC,
                type = "HC1",
                cluster = war) #stata robust 

stargazer(m1, m2, m3, m4,
          se = list(m1sr[,2], m2sr[,2],m3sr[,2], m4sr[,2]),
          covariate.labels = c("Corruption", "Initiation", "Target", 
                               "Liberal Democracy", "Military Executive Dependence",
                               "Capabilities", "Alliance Capabilities", "Troop quality",
                               "Terrain", "Strategy*Terrain","Strategy 1", "Strategy 2",
                               "Strategy 3", "Strategy 4"),
          dep.var.labels= c('Downes (2009) Win/Draw/Lose',
                            'Reiter and Stam (2002) Win/Lose',
                            'Downes (2009) Win/Lose',
                            "Cochran and Long (2017) Battle Performance"))




#####PREDICTED PROBABILITIES: SIMULATED on Model 3

Sims = 1000
newdata1 <-  expand.grid(constant = 1,
                         relative_corr = seq(min(LEDs_data_f$relative_corr), max(LEDs_data_f$relative_corr), 0.001), 
                         initally = 1,
                         targally = 0,
                         v2x_libdem = mean(LEDs_data_f$v2x_libdem),
                         v2x_ex_military = mean(LEDs_data_f$v2x_ex_military),
                         concap = mean(LEDs_data_f$concap),
                         capasst = mean(LEDs_data_f$capasst),
                         qualrat = mean(LEDs_data_f$qualrat),
                         terrain = mean(LEDs_data_f$terrain),
                         straterr = mean(LEDs_data_f$terrain)*3,
                         strat1 = 0, 
                         strat2 = 0, 
                         strat3 = 1, 
                         strat4 = 0)
                         

myCoefs = coef(m3)
myVCOV = vcov(m3)

# Now simulation the coefficients and their covariances
bhat = mvrnorm(Sims, myCoefs, myVCOV)

# bhat is 1000 by p t(xNew) is p by 2, calc ystar, linear index 
#  for the xNew scenarios
ystar = bhat %*% t(newdata1)  

# push ystar, the linear index, through the inverse logit function
# to create distribution of predicted probabilities
yhat = inv.logit(ystar)
# calc mean for the two scenarios
mean<-apply(yhat, 2, mean)
# calc quantile of interest 
quant<-apply(yhat, 2, quantile, probs=c(.05, .95))

preds<-t(rbind(mean,quant))

preds<-data.frame(cbind(newdata1,preds))

m3_plot <- ggplot(preds, aes(x=relative_corr, y=mean))+
  geom_ribbon(aes(ymin=`X5.`,ymax =`X95.`, fill=aa.index), alpha=.4, colour = NA,show.legend = F, fill=c("grey81")) +
  geom_line() +
  geom_rug(data=LEDs_data_f, aes(x = relative_corr,y = NULL)) + 
  labs(title="Model 3: Battle Performance", 
    y="Predicted Probability", x="Relative Corruption") +
  theme_bw()
#geom_rug(data.frame(full_data), mapping = aes(x=qmean1,y=1, group = NULL, color = NULL), show.legend = F, sides = "b")

m3_plot

#####PREDICTED PROBABILITIES: SIMULATED on Model 4

Sims = 1000
newdata1 <-  expand.grid(constant = 1,
                         relative_corr = seq(min(LEDs_data_f$relative_corr), max(LEDs_data_f$relative_corr), 0.001), 
                         initally = 1,
                         targally = 0,
                         v2x_libdem = mean(LEDs_data_f$v2x_libdem),
                         v2x_ex_military = mean(LEDs_data_f$v2x_ex_military),
                         concap = mean(LEDs_data_f$concap),
                         capasst = mean(LEDs_data_f$capasst),
                         qualrat = mean(LEDs_data_f$qualrat),
                         terrain = mean(LEDs_data_f$terrain),
                         straterr = mean(LEDs_data_f$terrain)*3,
                         strat1 = 0, 
                         strat2 = 0, 
                         strat3 = 1, 
                         strat4 = 0)


myCoefs = coef(m4)
myVCOV = vcov(m4)

# Now simulation the coefficients and their covariances
bhat = mvrnorm(Sims, myCoefs, myVCOV)

# bhat is 1000 by p t(xNew) is p by 2, calc ystar, linear index 
#  for the xNew scenarios
ystar = bhat %*% t(newdata1)  

# push ystar, the linear index, through the inverse logit function
# to create distribution of predicted probabilities
yhat = inv.logit(ystar)
# calc mean for the two scenarios
mean<-apply(yhat, 2, mean)
# calc quantile of interest 
quant<-apply(yhat, 2, quantile, probs=c(.05, .95))

preds<-t(rbind(mean,quant))

preds<-data.frame(cbind(newdata1,preds))

m4_plot <- ggplot(preds, aes(x=relative_corr, y=mean))+
  geom_ribbon(aes(ymin=`X5.`,ymax =`X95.`, fill=aa.index), alpha=.4, colour = NA,show.legend = F, fill=c("grey81")) +
  geom_rug(data=LEDs_data_f, aes(x = relative_corr,y = NULL)) + 
  geom_line() +
  labs(title="Model 4: Win/Lose", 
    y="Predicted Value", x="Relative Corruption") +
  theme_bw()
#geom_rug(data.frame(full_data), mapping = aes(x=qmean1,y=1, group = NULL, color = NULL), show.legend = F, sides = "b")

m4_plot 

###Combine Plots
grid.arrange(m3_plot,m4_plot,
             layout_matrix = rbind(c(1,2)
             ))


#######Look at Chechnya Battle Deaths#####
ucdp_ged <- read.csv(here::here("GEDEvent_v23_1.csv"))
ucdp_chech <- ucdp %>% filter((dyad_new_id == 852 & year <=1997))
unique(ucdp_chech$side_b)
unique(ucdp_chech$side_a)
sum(ucdp_chech$deaths_unknown)
sum(ucdp_chech$deaths_civilians)
sum(ucdp_chech$deaths_a)
sum(ucdp_chech$deaths_b)
sum(ucdp_chech$best)
sum(ucdp_chech$low)
sum(ucdp_chech$high)
